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ABSTRACT 

We carry out global three-dimensional radiation hydrodynamical simulations of self- 
gravitating accretion discs to determine if, and under what conditions, a disc may 
fragment to form giant planets. We explore the parameter space (in terms of the 
disc opacity, temperature and size) and include the effect of stellar irradiation. We 
find that the disc opacity plays a vital role in determining whether a disc fragments. 
Specifically, opacities that are smaller than interstellar Rosseland mean values promote 
fragmentation (even at small radii, R < 25AU) since low opacities allow a disc to cool 
quickly. This may occur if a disc has a low metallicity or if grain growth has occurred. 
With specific reference to the HR 8799 planetary system, given its star is metal-poor, 
our results suggest that the formation of its imaged planetary system could potentially 
have occurred by gravitational instability. We also find that the presence of stellar 
irradiation generally acts to inhibit fragmentation (since the discs can only cool to the 
temperature defined by stellar irradiation). However, fragmentation may occur if the 
irradiation is sufficiently weak that it allows the disc to attain a low Toomre stability 
parameter. 

Key words: accretion, accretion discs - planetary systems: formation - planetary 
systems: protoplanetary discs - gravitation - instabilities - hydrodynamics 



1 INTRODUCTION 

There are two ways in which giant planets have been hy- 



pothesised to form: core accretion ( Safronov|1969 Goldreich 


|& Ward||1973| |Pollack et al. 


|1996|) and gravitational insta- 


bility (Cameron 1978 Boss 


1997 review by Durisen et al. 



The former has been favoured but historically has 
had difficulties for two reasons: the first is a temporal is- 
sue where the timescales required to form planets may be 
too large such that the gas in the disc is depleted before 
the gas giant planet is fully formed. Secondly, simulations of 
this method model the growth of planets typically starting 
with kilometre-sized planetesimals, but while the growth of 
particles from small grains to metre-sized objects appears to 
be straight-forward, how to get from metre-sized objects to 
kilometre-sized planetesimals is unknown. 

Gravitational instabihty, on the other hand, eliminates 
the timescale problem, forming gas giant planets in < 0(10*) 
years. Such planets may not have solid cores, which may well 
be the case for Jupiter ( jSaumon fc Guillot|2004[ ), though the 
capture of solid material to form a core after the fragment- 
ing stage in a gravitationally unstable disc has also been 
proposed ( [Helled, Podolak fc Kovetz||2006[ ). However, since 
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gravitational instability is not thought to operate close to 
the central star, it has not been thought to be the dominant 
mechanism by which giant planets form as it was unable to 
describe the observations of close-in giant planets. Recent 



Marois et al. 


2008 1 encourage 


fc Bjorkman 


(|2009|l have arg 



haut b and at least the outer planet of the HR 8799 system 
could have formed by gravitational instability as the cooling 
timescales are likely to be small enough such that fragmen- 
tation is possible. Moreover, discs in their early stages are 
thought to be massive ( [Eisner fc Carpenter|20"06| suggesting 



that gravitational instability must play a role in the evolu- 
tion of a disc in the late Class I and early Class H stages. It 
has also been proposed that core accretion may be a method 
by which planets may form at small radii (~ O(IO)AU) 
whilst gravitational instability may be the mechanism by 
which planets may form at larger radii (> O(IQO)AU) (e.g. 
[Boley[[2009| ) though a hybrid scenario of forming gas giants 
in the same system by both core accretion and gravitational 
instability has yet to be modelled. 

There are two quantities that have historically been 
used to determine whether a disc is likely to fragment. The 
first is the stability parameter ( Toomre|1964 |, 
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Q = 



ttEG' 



(1) 



where Cs is the sound speed in the disc, Kop is the epicyclic 
frequency, which for Keplerian discs is approximately equal 
to the angular frequency, Q, E is the surface mass density 
and G is the gravitational constant. Therefore, once the sur- 
face mass density and the rotation of the disc have been 
established, the stability is purely dependent on the disc 
temperature. Toomre (19641 showed that for an infinitesi- 



mally thin disc to fragment, the stability parameter must 
be less than a critical value, Qcrit ~ 1- 



Gammie (2001 1 showed that in addition to the stability 



criterion above, the disc must cool at a fast enough rate. Us- 
ing shearing sheet simulations, he showed that if the cooling 
timescale can be parametrized as 



where 



^cool 



/ dltcool 



(2) 



(3) 



u is the internal energy and ducooi/dt is the total cooling 
rate, then for fragmentation we require P < 3, for a ratio 
of specific heats 7 = 2 (in two dimensions). Rice, Lodato 



& Armitage ( 2005 1 carried out three-dimensional simula- 



tions using a smoothed particle hydrodynamics (SPH) code 
and showed that this cooling parameter is dependent on the 
equation of state so that fragmentation can occur for discs 
with 7 = 5/3 and 7/5 if /3 < 7 and 13, respectively. [Clarke, | 
Harper-Clark fc Lodato| ( |2007 1 showed that the critical value 



of P (below which fragmentation will occur if the stability 
criterion is met) may depend on the disc's thermal history: 
if the timescale on which the disc's cooling timescale is de- 
creased is slower than the cooling timescale itself (i.e. a grad- 
ual decrease in /3) then the critical value may decrease by 
up to a factor of 2. More recently, |Cossins, Lodato fc Clarke] 
( 2010 1 showed that the critical value varies with the temper- 



ature dependence of the cooling law. In general though, the 
critical value is thought to be of the order of the dynamical 
timescale. 

The above fragmentation criteria are based on the as- 
sumption that the dominant form of dissipation in the disc 
is due to internal heating processes. Previous simulations 
without external irradiation have considered isolated discs 



with simple cooling prescriptions (e.g. Lodato & Rice 



and with radiative transfer (e.g. Boss|20(5T 



Cai et al 



20041 



2006 



[Mayer et a"L]|2007| . [Johnson fc Gammie| ( [2003[ ) suggested 
that discs with external irradiation are likely to be effectively 
isothermal and can therefore be treated as such. IMatznerl 
[fc Levin] ( |2005| l analytically considered externally irradiated 
discs and concluded that stellar irradiation quenches frag- 
mentation. Cai et al. ( 2008 1 carried out simulations with 



external irradiation and found that their discs are more re- 
sistant to fragmentation and proposed that these results may 
be extended to discs with stellar irradiation. IStamatellos fcl 
[Whitworth[ ( [2008[ ) also carried out simulations taking into 
account the effects of stellar irradiation and found this to 



be a stabilising factor. Rafikov (20091 analytically explored 



fragmentation in gravitoturbulent discs including the effects 
of stellar irradiation and suggested that fragmentation can 
only occur beyond ~ 120AU. Dodson-Robinson et al. (20091 



carried out a linear stability analysis on irradiated discs to 
show that gravitational instability takes place for systems 



with a large disc to star mass ratio. However, whilst frag- 
mentation in gravitationally unstable discs may be less likely 
than previously thought, we still do not know in what sit- 
uations discs may fragment when modelling them realisti- 
cally with radiative transfer and by considering the effects of 
stellar irradiation. It is therefore important to deduce when 
fragmentation may occur when simulating discs with more 
detailed energetic conditions, and just how realistic or un- 
realistic fragmentation is in real discs. 

[Boss[ ( [2002[ ) carried out simulations of gravitationally 
unstable discs and varied the opacity from O.lx to 10 x the 
Rosseland mean opacities and found that the fragmentation 
results were insensitive to the dust grain opacity. However, 
given that a reduced opacity is more likely to allow energy 
to stream out of a disc more easily causing it to cool and 
promote fragmentation, whilst in a high opacity disc the 
converse is true, it is interesting to consider what opacity 
values allow and do not allow fragmentation. Given that a 
disc's opacity gives somewhat an indication of how metal- 
rich it is or how large or small the grain sizes are, we may 
then make preliminary conclusions on the disc conditions 
that are likely to promote fragmentation, which is a key 
focus of this paper. 

In this paper, we model the evolution of massive self- 
gravitating discs using a global three-dimensional SPH code 
including radiative transfer and the effects of stellar irradia- 
tion. In particular, we explore the parameter space in terms 
of dust opacity, disc temperature and size in order to scope 
out if, and under what conditions, a self-gravitating disc may 
fragment. In Section [2j we describe the code used to carry 
out our simulations. In Section[3j we outline our simulations, 
including the disc setup and discuss the parameter space. In 
Section [4] we present our results, while we compare with 
previous studies and make conclusions in Sections [5] and [6] 
respectively. 



2 NUMERICAL SETUP 

Our simulations are carried out using an SPH code origi- 



nally developed by Benz ( 1990 1 and further developed by 



Bate, Bo nnell fc Price[ ([1995[), [Whitehouse, Bate fc Mon 
aghan (2005[), [Whitehouse fc lBate[ ( [2006[ ) and [Price fc Bate 
(2007). It is a Lagrangian hydrodynamics code, ideal for 
simulations that require a large range of densities to be fol- 
lowed, such as fragmentation scenarios. The version used 
to carry out the simulations presented here includes radia- 
tive transfer using the flux-limited diffusion approximation 
( Whitehouse et al.|2005| [Whitehouse fc Bate|2006[ ) with two 
temperatures: that of the gas and that of the radiation field. 

In order to model shocks, SPH requires artificial viscos- 
ity. We use a common form of artificial viscosity by |Mon-| 
aghan & Gingold (1983), which uses the parameters asPH 



and PspB.- A corollary of including artificial viscosity is that 
it adds shear viscosity and causes dissipation. If this vis- 
cosity is too large, the evolution of the disc may be driven 
artificially, while if it is too small, it will lead to inaccu- 
rate modelling of shocks (jBatej^l995|) . We have tested vari- 
ous values of the SPH parameters and find that a value of 
aspH ~ 0.1 provides a good compromise between these fac- 
tors. Since typically, Psph ~ 2aspH, we choose aspH and 
PsPH to be 0.1 and 0.2, respectively. These are the same val- 
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ues as those implemented by Lodato & Rice (20041. How- 



ever, we have chosen to fix our artificial viscosity parame- 
ters whereas Lodato & Rice' ( 2004 ) used the 'Balsara (19951 



switch which reduces the artificial viscosity in places where 
the ratio of the divergence to the curl of the velocity field 
is small i.e. in regions of strong vorticity. The reasoning be- 
hind our choice is two-fold: (i) the Balsara switch reduces 
the viscosity in shearing flows but given that we have both 
shocks and shearing flows, the reduction does not handle the 
shocks well; (ii) we want to ensure a controlled test which is 
not possible if the numerical viscosity is constantly chang- 
ing. 



atmosphere temperature is not set by the disc itself but by 
some other process such as stellar or external irradiation. 
The particles forming the disc boundary are variable since 
some particles may move across the boundary from the op- 
tically thin to the optically thick region and vice versa. We 
have carried out tests to see what the effects of the exact 
location of the vertical boundary has on the disc energet- 
ics and have found that provided the boundary is located 
in the optically thin surface region of the disc and provided 
that a layer of boundary particles exists over the entire disc 
surface, the disc energetics remain unchanged. 



2.1 Opacity &: equation of state 

The discs simulated are assumed to be in local thermal equi- 
librium and the opacities are assumed to be grey Rosseland 
mean values. The opacities are based on the opacity tables 
of interstellar molecular dust grains pro duced by [Pollack, | 
McKay fc Christofferson| (|1985[), and on |Alexander| (|1975[) 



for the higher temperature gaseous contributions. A sum- 
mary of how the opacity tables were used is available in 



Whitehouse & Bate ( 2006 1 . For non-interstellar opacity sim- 



ulations, we scale these values by the required factor (see 
Table [l] for simulation details). The details of how this was 
done can be found in Section 2.4 of'Ayhffe fc Bate| (["2009 1 . 



The equation of state used in these simulations assumes 
that the gas is composed of hydrogen (70%) and helium 
(28%) and includes dissociation of molecular hydrogen, ion- 
isation of both hydrogen and helium, and the rotational 
and vibrational modes of molecular hydrogen, as done so 



by Whitehouse & Bate (20061, though the equation of state 



has been corrected following Boley et al. ( 2007 1. It omits the 
contribution due to metals. Our discs are cold enough such 
that the rotational and vibrational degrees of freedom are 
not excited so that effectively, 7 = 5/3. Our models assume 
that the presence of dust only affects the value of the opacity 
in the disc. We do not include the effects of dust on the disc 
dynamics, nor do we consider its effects on the equation of 
state. 

2.2 Stellar irradiation 

As mentioned earlier, radiative transfer is simulated using 
the flux-limited diffusion approximation. We use a two-layer 
approach to simulate the midplane and surface regions. In 
order to model the energy loss from the disc surface, a 
boundary layer of particles maintains a fixed temperature 
profile such that any energy that is passed to these boundary 
particles is effectively radiated away. The vertical location of 
the boundary between the optically thick part of the disc and 
the optically thin atmosphere is located at the maximum 
of: the height above the midplane where the optical depth, 
r = 1, or the height above the midplane where Zh — 1.75//, 
where H is the isothermal scale height given hy H — Cs/Q. 
Therefore, the boundary is at the vertical position where 
the optical depth, t < 1, and our choice ensures that the 
"bulk" of the disc (i.e. consisting of particles that lie closer 
to the disc midplane) is simulated using radiative transfer, 
rather than simplified energetic calculations. The remaining 
"boundary" particles have their temperature held at a con- 
stant temperature profile. Physically, we assume that the 



3 SIMULATIONS 

Table [l] summarises the parameters and fragmentation re- 
sults of the simulations presented here. Each simulation was 
run either beyond the point at which the disc attained a 
steady state (for > 6 outer rotation periods, ORPs), or un- 
til it fragmented (defined as regions which are at least three 
orders of magnitude denser than their surroundings). The 
magnitude and profile of the boundary temperatures are the 
same as those of the initial discs so that the discs start in 
thermal equilibrium with their vertical boundaries (atmo- 
spheres). We expect the discs to heat up initially due to 
work done and viscous heating. Following an initial tran- 
sient phase, the bulk of the disc may or may not cool to 
re-establish thermal equilibrium with the boundary. We dis- 
cuss a Reference case first before turning our attention to 
exploring the parameter space. 



3.1 Reference case 

Our reference disc is set up in exactly the same way as 
a disc simulated by Lodato & Rice (20041: a 1 Mq star 
with a 0.1 Mq disc made of 250,000 SPH particles, span- 



ning 0.25 5S -R 5S 25AU. The initial surface mass density 
and temperature profiles of the disc are E oc and 
T oc R~ 3 , respectively. The magnitudes of these are set 
such that the Toomre stability parameter (equation [T| at 
the outer edge of the disc, Qmin = 2. This gives an aspect 
ratio, H/R ~ 0.05. We model the 1 Mq star in the centre of 
the disc using a sink particle ( Bate et al.|1995 I . At the inner 
disc boundary, particles are accreted onto the star if they 
move within a radius of 0.025 AU of the star or if they move 
into 0.025 ^ R < 0.25AU and are gravitationally bound to 
the star. At the outer edge, the disc is free to expand. 

The initial Reference disc is a Toomre stable disc. Given 
that the boundary temperature profile is the same as that 
of the initial disc and that Qmin = 2, we do not expect this 
disc to fragment. However, we use the Reference disc as a 
fiducial case. In particular, we are concerned with the cool- 
ing rates that are present in the disc once it is in an equi- 
librium state. We emphasize our use of terminology here: 
when referring to the disc being in thermal equilibrium with 
the boundary, we are describing the bulk of the disc being a 
similar temperature to the disc boundary (which is assumed 
to be determined by stellar irradiation), whereas an equi- 
librium state disc refers to the dissipative and cooling rates 
being balanced such that the Toomre stability profiles do 
not change with time. 
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Simulation 


Disc radius 


Opacity 


Qmin 


Fragment? 


Fragmentation 


name 


[AU] 


scaling factor 






time 


Reference 


25 


1 


2 


n 


- 


KappalO 


25 


10 


2 


n 


- 


KappaO.l 


25 


0.1 


2 


n 




KappaO.Ol 


25 


0.01 


2 


n 




Qminl 


25 


1 


1 


n 




QminO.75 


25 


1 


0.75 


n 




QminO. 75-KappaO . 1 


25 


0.1 


0.75 


n 




Qmin0.75-Kappa0.01 


25 


0.01 


0.75 


y 


9.7 ORPs 


QminO. 5 


25 


1 


0.5 


y 


1.8 ORPs 


L-Qminl 


300 


1 


1 


n 




L-Qminl-KappalO 


300 


10 


1 


n 




L-Qminl-KappaO. 1 


300 


0.1 


1 


y 


1.5 ORPs 


L-QminO.75 


300 


1 


0.75 


y 


< 1 ORPs 



Table 1. Summary of the simulations carried out. The opacity scalings refer to multiples of interstellar Rosseland mean opacity values 
as described in Section [2] Qmin refers to the minimum value of the Toomre parameter (at the outer edge of the disc) at the start of the 
simulation. 



3.2 Exploring the parameter space 

Given that a motivation of this work is to determine if, 
and under what circumstances, fragmentation in realisti- 
cally modelled self-gravitating discs may occur, we explore 
the parameter space in a number of ways. One parameter is 
the opacity: we re-run the Reference simulation with opacity 
values scaled to 10 x, O.lx and 0.01 x the interstellar opac- 
ity values (simulations KappalO, KappaO.l and KappaO.Ol, 
respectively) . This may be equivalent to a disc with differing 
metallicities or grain sizes. We make two assumptions here: 
i) the change in metallicity does not affect the equation of 
state (as described in Section 2.1 1 and ii) there are no spa- 
tial or temporal variations in the grain size distributions. As 
with the Reference case, these discs are simulated purely to 
analyse the energetics since we do not expect these discs to 
fragment. 

We then choose to explore the initial and boundary tem- 
perature conditions by decreasing the magnitude of the disc 
temperature whilst maintaining the same surface mass den- 
sity as the Reference case. We do this by changing the initial 
Toomre stability parameter profiles such that Qmin = 1, 0.75 
and 0.5 (simulations Qminl, QminO. 75 and QminO. 5, respec- 
tively). This is equivalent to reducing the disc aspect ratios 
to H/R ~ 2.2 X 10~^ 1.7 X 10"^ and 1.1 x 10"^ respec- 
tively. We reiterate that the boundary temperature is the 
same as the temperature of the initial disc and hence this 
setup not only changes the disc temperature profile, but it 
also changes the boundary temperature profile. 

Furthermore, we consider a combination of the above 
factors by simulating discs with Qmin = 0.75 and opacities 
that are O.lx and 0.01 x the interstellar opacity values. 

The unfavourable conditions for fragmentation at small 
radii have been discussed at great length in the past (e.g. 
[2005| JStamatellos fc Whitworth|[2008l |Boley|[2509l 



Rafikov 



Rafikov 



2009 



Clarke 



2009 1 . We therefore expand our param- 



eter space to include discs that are a factor of 12 larger with 
a radii range of 3 ^ 7? ^ 300AU. These discs have the same 
mass as the 25 AU discs and are set up so that Qmin = 1. We 
simulate three different opacity values (lx,10x and O.lx 
the interstellar Rosseland mean opacities). In addition, we 



also simulate a large disc with Qmin = 0.75 with interstellar 
opacity values. 

In order to keep these disc masses and initial Toomre 
stability profiles the same as the smaller 25AU discs, we 
require both the surface mass density and absolute tem- 
perature to be reduced. These discs are therefore not only 
larger, but also colder than their equivalent (in terms of ini- 
tial Toomre stability profiles) small discs. 



4 RESULTS 

The simulations have been analysed in three main ways: 

(i) we compare the azimuthally averaged Toomre stabil- 
ity profiles of the initial and final (or in the case of fragment- 
ing discs, shortly before fragmentation) discs which indicates 
whether the bulk of the discs were able to reach a state of 
thermal equilibrium with their boundary. The surface mass 
density does not change significantly throughout the simula- 
tions and hence changes in the Toomre stability parameter 
are due to changes in the disc temperature. This enables us 
to determine which discs are more likely to fragment. Note 
that we assume Kcp = in equation [l] 

(ii) we examine the timescale on which the discs cool 
(by considering the energy passed from the gas to the radi- 
ation within the disc as well as that which is assumed to be 
instantly radiated away from the disc surface by the bound- 
ary particles). In past simulations that have neglected the 
heating effects of stellar irradiation (e.g. Gammie|2001 Rice 
|et al.„2005[ ), the cooling, C, in a steady-state disc balances 
the heating due to gravitational stresses, Hqi, and the heat- 
ing due to artificial viscosity, i/^, such that 



(4) 



If the artificial viscosity is low, C ~ Hai- In this case, the 
cooling timescale in units of the orbital timescale, /3 (sec- 
tion [l|, can be related to the gravitational stress in the disc, 
aai ( |Gammie||2001[ ) : 

4 11 



Gammie (20011 and Rice et al. (2005) have shown that 
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the maximum gravitational stress that a disc can support 
is aci = 0.06, beyond which fragmentation will occur. In 
discs that do not take into account heating due to exter- 
nal irradiation, this condition is equivalent to requiring the 
cooling timescale in terms of the orbital timescale, /3, to be 
smaller than the critical values, described in section [l] for 
fragmentation. 

In our steady-state discs, not only does the cooling have 
to balance the heating due to the gravitational instabilities 
and the numerical viscosity, but it also has to balance the 
heating due to stellar irradiation, Hsi, such that: 



10 



C = Hgi+H^+ Hsi 



(6) 



In what follows, we calculate the parameter, tp, which we 
define to be the timescale on which the disc cools in units 
of the orbital timescale. Without irradiation, ip = j3. When 
including heating due to stellar irradiation, the tp parameter 
does not specifically tell us about the fraction of the cool- 
ing that balances the gravitational instabilities and therefore 
cannot be used directly to decide whether a disc should frag- 
ment or not. However, it does still give an indication as to 
what the cooling rate is in the discs which has been shown 
to be important when determining whether a disc is likely 
to fragment or not ( |Gammie|200T| |Rice et al.||2005[ ). 

(iii) we examine images of the discs for signs of fragmen- 
tation which we define as clumps in the disc which are at 
least three orders of magnitude denser than their surround- 
ings. 

4.1 Reference case 

Figure [T] shows the Toomre stability profile of the initial 
and final Reference disc. The disc is not able to cool rapidly 
enough in response to the internal heating and consequently 
the disc midplane becomes hotter than the boundary (top 
panel of Figure [2|. Figure [l] also shows the final Toomre sta- 
bility profile for the equivalent disc simulated by |Lodato fc| 
Rice ( 2004 \ which used simplified cooling rather than radia- 



tive cooling. We have included this comparison to emphasize 
how significant the differences can be between discs simu- 
lated using simplified cooling parameters and discs modelled 
not only with more detailed radiative cooling but also incor- 
porating the effects of stellar irradiation. 

4.2 Opacity effects 

Figure |3] shows the effect on the Toomre stability parame- 
ter of changing the opacity in the discs to 10 x, O.lx and 
0.01 X the interstellar Rosseland mean values (simulations 
KappalO, KappaO.l and KappaO.Ol, respectively). We see 
that decreasing the disc's opacity enhances its ability to 
reach a state of thermal equilibrium with the boundary 
since the radiation leaves the disc far more effectively and 
hence the disc cools faster. This is particularly evident in 
the cross-sectional plots showing the temperature structure 
in the vertical direction (Figure |2|. Figure [l] shows the ip 
profile of these simulations. This figure shows that a low 
opacity disc has a greater ability to radiate away the disc's 
energy. The decrease in the cooling timescale, tp, does not 
continue at very low opacities because the stellar irradiation 
sets the boundary temperature and therefore, the minimum 
disc temperature. 
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Figure 1. Azimuthally averaged values of the Toomre parameter 
at the start (solid line) and at a time t = 6.4 ORPs (dotted line) 
for the Reference simulation. The disc is unable to cool rapidly 
due to the internal heating and hence its end state is more stable 
than the initial disc. Also shown is the equivalent disc simulated 
by [Lodato &: Rice] | |2004[ short dashed line) which cools using 
simplified cooling (with /3 = 7.5) rather than by radiative cooling, 
and also does not consider the effects of stellar irradiation. The 
critical value of Qcrit = 1 is also marked. 




Figure 2. Logarithm of the gas temperature (in K) rendered in 
cross-sectional views of the Reference (top panel) and KappaO.Ol 
(bottom panel) discs at a time t = 6.4 ORPs. The surface of 
the Reference disc is clearly colder than the midplane whilst the 
midplane of the KappaO.Ol disc is closer to a state of thermal 
equilibrium with the boundary. Axis units are in AU. 



Figures |3] and |4] particularly show that the low opacity 
discs are able to cool fast enough to reach a state of thermal 
equilibrium with their boundaries. Thus, if the conditions 
were right (i.e. the boundary temperature was lower so that 
the Toomre stability parameter was able to reach Q < 1), 
the low opacity disc may fragment. We therefore turn our 
attention to the disc absolute temperature. 
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Figure 3. Azimuthally averaged values of the Toomre stability 
parameter at the start (sohd line) and at a time t = 6.4 ORPs, 
for the Reference (dotted line), KappalO (short dashed line), 
KappaO.l (long dashed line) and KappaO.Ol (dot-dashed line) 
simulations. The critical value of Qcrit = 1 is also marked. De- 
creasing the opacity causes the disc to cool more efficiently. 



Figure 5. Initial (thin lines) and final (at t = 6.4 ORPs; heavy 
lines) Toomre stability profiles for the Reference (solid line), 
Qminl (dotted line) and QminO.75 (dashed line) simulations. The 
critical value of Qcrit = 1 is also marked. None of the three discs 
which have interstellar opacity values can cool rapidly enough to 
maintain thermal equilibrium with their boundaries. 
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Figure 6. Surface density rendered image of the QminO.75 disc 
at a time t = 6.4 ORPs. At the end of the simulation, the disc has 
not fragmented despite initially being in a critical state because 
with interstellar opacities, it heats up. 



Figure 4. Cooling timescale, -0, profiles for the Reference (dotted 
line), KappalO (short dashed line), KappaO.l (long dashed line) 
and KappaO.Ol (dot-dashed line) simulations at a time t = 6.4 
ORPs. Decreasing the opacity causes the cooling rate in the disc 
to increase (and ip to decrease), but only until the disc can reach 
thermal equilibrium with its boundary. 



4.3 Colder discs 

Figure [5] shows the initial and final Toomre stability profiles 
for the Reference, Qminl and QminO.75 simulations. The 
interesting aspect about the QminO.75 and Qminl cases are 
that though our initial and boundary conditions are either 
unstable or marginally stable, the discs still do not fragment 
because they are unable to cool rapidly and heat up so that 
they end up being at or above the marginal state (e.g. Fig- 
ure |6|. The cooling timescales, ip, for the discs are as low as 
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Figure 7. Cooling timescale, tit, profiles for the Reference (solid 
line), Qminl (dotted line) and QminO.75 (dashed line) simulations 
at time t = 6.4 ORPs. With interstellar opacities, these discs arc 
simply not able to cool rapidly enough to obtain a ip value that 
is low enough for fragmentation. 



~ 20 (Figure [7| suggesting that for the discs to fragment, 
an efficient energy removing mechanism is needed such that 
the coohng timescale is lower than the QminO.75 (dashed 
hne) curve in Figure [t] 

For completeness, we also simulate a disc with 
Qmin = 0.5 and find that though this disc also heats up, it 
cannot heat fast enough to become Toomre stable before it 
fragments (due to its initial conditions). 



4.4 Low temperature, low opacity discs 

Until now, we have identified under what conditions discs are 
more likely to fragment (i.e. sub-interstellar opacities and 
cooler temperatures). Combining these conditions, we sim- 
ulate two further discs with outer Toomre stability parame- 
ters of 0.75 with O.lx and O.Olx interstellar opacity values 
(simulations Qmin0.75-Kappa0.1 and Qmin0.75-Kappa0.01, 
respectively). With such low opacity values, the QminO.75- 
KappaO.Ol disc would be equivalent to a low metallicity disc 
or a disc with grain sizes ranging between millimetre and 
centimetre sizes which is realistic given current observations 



( Calvet et al.|2002| |Testi et al.|[2003} [Rodmann et al.|[2006 



Lommen et aL||2007 l 

We can see from Figure |8] that the disc with O.lx inter- 
stellar opacity values is slightly cooler than the case with in- 
terstellar opacities, but is still unable to cool rapidly enough 
to maintain a state of thermal equilibrium with its bound- 
ary. However, the disc with O.Olx interstellar opacity values 
is able to since its Toomre stability profile at the end of the 
simulation is close to its initial value. We find that the disc 
in the latter simulation does indeed fragment (Figure [9|, 
though it takes ~ 9.7 ORPs to do so. This is because: 

(i) the disc goes through a transient phase where it heats 




10 15 
Radius [AU] 

Figure 8. Toomre stability profiles at t = 6.4 ORPs for the 
QminO.75 (short dashed line) and Qmin0.75-Kappa0.1 (long 
dashed line) simulations and at t = 8 ORPs (just before it begins 
to fragment) for the Qmin0.75-Kappa0.01 (dotted line) simulation 
in comparison to the initial (solid line) Toomre stability profile 
for these discs. The critical value of Qcrit = 1 is also marked. 
The lower opacity disc cools rapidly enough to attain a state of 
thermal equilibrium with its boundary. It eventually fragments at 
t = 9.7 ORPs 




Figure 9. Surface density rendered image of the fragmented 
Qmin0.75-Kappa0.01 disc at a time of t = 10.5 ORPs. The disc 
not only requires reduced irradiation, but also low opacities are 
essential for it to cool rapidly enough to fragment. 



up since the initial disc is not quite in hydrostatic equilib- 
rium. 

(ii) although the radiation is able to leave quickly, the 
disc is enveloped in a thermal blanket due to the effects of 
stellar irradiation, thus causing the disc to cool more slowly. 

(iii) as the disc midplane cools, its cooling rate also 
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Figure 10. Cooling timescale, ip, profiles for the QminO.75- 
KappaO.Ol simulation (dotted line) at t = 8 ORPs, just before it 
begins to fragment, in comparison to the QminO.75 (short dashed 
line), KappaO.Ol (long dashed line) and Qmin0.75-Kappa0.1 (dot- 
dashed line) simulations at t = 6.4 ORPs. The fragmenting disc 
has a low value (Ri 15) in the outer part of the disc where it 
fragments. 



Figure 11. Toomre stability profiles for simulations L-Qminl 
(dotted line), L-Qminl-KappalO (short dashed line) and L- 
Qminl-KappaO.l (long dashed line) at a time of t = 6.4 ORPs 
as well as the boundary profile for these discs (solid line). The 
critical value of Qcrit = 1 is also marked. Larger discs do not re- 
quire as low opacities as smaller discs to attain a state of thermal 
equilibrium with their boundaries. 



decreases since the temperature gradient between the disc 
midplane and surface becomes smaller. 

Figure [TO] shows the cooling timescale, i/i, of simulation 
Qmin0.75-Kappa0.01 (1.7 ORPs prior to fragmentation) in 
comparison to simulations KappaO.Ol (which has Qmin = 2), 
QminO.75 (which uses interstellar opacities) and QminO.75- 
KappaO.l. With a reduced opacity, the Qmin0.75-Kappa0.01 
disc cools on a timescale fast enough such that it fragments 
i.e. it cools fast enough so that the heating due to gravi- 
tational instabilities, numerical viscosity and stellar irradi- 
ation do not heat the disc significantly above the boundary 
temperature. This figure shows that the cooling timescale 
-i/j « 15 in the outer parts of the disc, where fragmentation 
occurs. We emphasise the need to express caution when con- 
sidering the absolute value of i/; in a non-steady state frag- 
menting disc, since it is dependent on the disc temperature 
which is constantly changing. We discuss this in detail in 
Section [5] A key interesting aspect about this simulation 
is that the fragments form even though the disc is small 
(R < 25 AU), a result contrary to what has been suggested 
in the past (see Section [5] for further discussion). 




Figure 12. Surface density rendered image of the large, low opac- 
ity disc, L-Qminl-KappaO.l, at a time of t = 2.9 ORPs. A low 
opacity is also required for a largo disc to fragment, though the 
opacity does not have to be as low as the 25 AU discs. 



4.5 300 AU discs 

We now turn our attention to larger discs that are usually 
considered more likely to promote fragmentation. 

Simulations L-Qminl-KappalO, L-Qminl and L- 
Qminl-KappaO.l show that though it is easier for larger 
discs to maintain a state of thermal equilibrium with their 



boundaries (Figure 111, higher opacity discs still struggle to 



do so, thus implying that low opacities are still important 
for discs to fragment. 

However, despite simulating a large disc with a marginal 
Toomre stability profile (Qmin ~ 1; simulation L-Qminl), 
we are still unable to see any fragmentation for discs with 
interstellar opacities unless we simulate discs whose frag- 
mentation is determined by the initial conditions (e.g. Sim- 
ulation L-QminO.75 which fragments in t < 1 ORP). We 
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find that though the opacities do not have to be as small as 
those in the 25 AU discs (0.01 x interstellar values), never- 
theless, we still require discs with sub-interstellar opacities 
(O.lx interstellar values). Figure 12 shows the disc of simu- 



tions A1-A5 of Clarke ( 2009 1 show that if the opacity is 



lation L-Qminl-KappaO.l which fragmented after more than 
an outer rotation period of evolving. Though for fragmen- 
tation to occur, the opacity in this disc does not have to 
be decreased as much as for the smaller disc simulated in 
Qmin0.75-Kappa0.01, since it is colder, the grain sizes still 
correspond to sizes ranging between millimetre and centime- 
tre which, as mentioned earlier, is entirely realistic. 



5 COMPARISON WITH PREVIOUS WORK 

Our simulations show that contrary to past studies, it is 
possible for discs to fragment at small radii (< 25 AU) if the 
disc temperature and the opacity are low enough. The latter 
point may be the case if the disc is metal-poor or if grain 
growth has occurred to produce grains that are larger than 
interstellar sizes. The larger grain size is certainly a reason- 
able assumption since there is evidence for grains of up to 
centimetre sizes in discs ( Calvet et al.|20d2 Testi et al. 2003 



Rodmann et al.|2006 Lommen et al. 2007 1, which could pro- 
vide the low opacities necessary for our fragmenting discs. 
With respect to disc metallicity, it is well known that core 



accretion is not as efficient in metal-poor environments ( Ko- 
rnet et al.|2005 I. However, we find that gravitational insta- 



bility is enhanced in such conditions. Boss (20021 carried 
out simulations of gravitationally unstable discs with vary- 
ing opacities (from O.lx to 10 x the interstellar Rosseland 
mean values) in order to explore fragmentation in different 
metallicity discs and found that the fragmentation results 
were insensitive to the dust grain opacity. He also found that 
the disc midplanes could radiate energy to the disc surfaces 
very rapidly as the timescale for temperature equilibration 
by radiative diffusion between the disc midplane and the sur- 
face was smaller than the orbital timescale regardless of the 
opacity scaling. However, for our discs the converse is true 
and we find that the opacity, and hence metallicity, does 
play a part in the likelihood of fragmentation. Our results 
suggest that gravitational instability may be the dominant 
giant planet formation mechanism in low metallicity discs. 

Our small fragmenting discs are also in contrast to pre- 
vious work which has suggested via simulations (e.g. |Sta- 



matellos & Whitworth| 


2008| |Boley| |2009|) and analytical 


work (e.g. Rafikov 2005 


2009 Clarke 2009 1 that fragmen- 



tation at such small radii is not possible. However, this is 
due to the lower opacities in our discs. I Rafikov] ( [20091 ) ana- 
lytically explored self-gravitating discs including the effects 
of stellar irradiation and suggested that fragmentation in- 
side of ~ 120AU is not possible. However, firstly he uses 
interstellar Rosseland mean opacities (which is not an un- 
reasonable assumption to initially make) and secondly he 
assumes that the fragmentation boundary occurs when the 
"effective a parameter" (of the form of Shakura & Syun 



yaev|1973 1 due to the gravitational torques oqi ~ 1- Other 
authors (e.g. |Clarke|2009] ) assume the fragmentation bound- 
ary occurs when ogi ~ 0.06 (Gammie 2001 Rice et al.|2005 1 
which occurs at a radius R ~ 70 AU. This is still at a larger 
radii than our small fragmenting disc, though [Clarke| ( [2b09[ ) 
also assumes interstellar Rosseland mean opacities. Equa- 



decreased by two orders of magnitude, as is the case for our 
small fragmenting disc, and using aai ~ 0.06 as the value at 
which fragmentation can occur, the radius outside of which 
fragmentation occurs is at i? ~ 24 AU. This is much smaller 
than the canonical value of ~ 70 AU, although still not quite 
as small as i? « 15 AU at which we find fragmentation. 
Given that the work of Clarke ([2009!) was analytical and 
while we have performed direct global three-dimensional ra- 
diation hydrodynamical simulations, this level of agreement 
is reasonable. 

We can also compare our results to the analytical work 



of Rafikov (20051. Based on a combination of the stability 



and cooling criteria of Toomre (1964) and Gammie (2001), 



respectively, Rafikov ( 2005 I analytically derived a constraint 
on the disc temperature that is required for planet formation 
via gravitational instability. Using this constraint (equation 
5 of Rafikov 2005 1 and taking into account the reduced opac- 



ity in our discs, the Rafikov fragmentation conditions for the 
surface mass density and temperature in the disc become 



QcritTrG 
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and 



(7) 



(8) 



respectively, where k is the Boltzmann constant, /i is the 
mean particle mass, a is the Stefan-Boltzmann constant, Ko 
is the constant of proportionality for the Rosseland mean 
opacity expression for this temperature and low opacity 
regime given hy k — KoT^ and has a value Ko = 2 x 
10~®g cm~^ K~^, and ^ = 2/3(7 ~ 1) s-nd Eilso absorbs any 
0(1) factors that have not been accurately considered in a 
proper calculation of the cooling time. Rafikov ( 2005 1 as- 
sumes that ~ 1. We find that our fragmenting disc is in 
agreement with these conditions to with a factor of ~ 2. As 



with the Clarke ( 2009 ) comparison, this level of agreement 
is reasonable given that our simulations are global, three- 
dimensional and use a realistic radiative transfer method. 
Consequently, our simulation results are consistent with pre- 
vious analytical work, but they emphasise the importance of 
opacity in determining the radius outside of which fragmen- 
tation may occur. 

Another difference between our simulations and those 
in the past that have used simplified cooling is that the heat- 
ing in previous simulations has been dominated by internal 
heating processes i.e. heating due to gravitational instabil- 
ity and viscous processes, whereas our simulations involve 
additional external heating due to stellar irradiation. This 
additional thermal blanket causes the Toomre stability pa- 
rameter to remain high, thus inhibiting fragmentation (con- 
sistent with the results of ICai et 8171120081 and IStamatellosI 
[fc Whitworth|2008[ ). However, in scenarios where the stellar 
irradiation is not so strong such that the Toomre stability 
parameter is small, fragmentation is possible, but only if the 
opacity is also decreased such that the disc can cool easily. 

The cooling timescale, i/i, parameter cannot give any 
indication as to what the gravitational stress in a disc is 
since it incorporates the cooling in response to the heating 
due to stellar irradiation as well as the gravitational stresses. 
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However, it is still interesting to compare the values in our 
discs to the critical values of P for fragmentation obtained 
from simulations using more simplified energetics. Previous 
simulations using simplified cooling have suggested that for 
discs with 7 — 5/3, fragmentation requires that the total 
cooling timescale in units of the orbital timescale, /3 < 7 



(Rice et al. 20051. Our non-fragmenting discs which have 



reached a steady state are consistent with this result since 
for these discs -0 > 7. For our fragmenting disc, we find 
that the cooling timescale in units of the orbital timescale 
can be as much as ~ 15. However, we express caution 
when interpreting this result. It is important to note that the 



Gammie ( 2001 1 and Rice et al. ( 2005 I conditions indicate the 



minimum cooling time (and hence a maximum gravitational 
stress) that a disc can support without fragmenting. In our 
fragmenting disc, the temperature continues to change as 
the disc cools, and hence the cooling rate is not constant 
since it is temperature dependent. 



disc's temperature does not increase significantly (due to in- 
ternal dissipation) above the boundary temperature. 



7 ACKNOWLEDGMENTS 

We thank the anonymous referee whose comments greatly 
improved the clarity of the paper. The calculations reported 
here were performed using the University of Exeter's SGI 
Altix ICE 8200 supercomputer. The disc images were pro- 
duced using SPLASH ( |Price|2007| . MRB is grateful for the 
support of a EURYI Award which also funded FM. This 
work, conducted as part of the award 'The formation of 
stars and planets: Radiation hydrodynamical and magne- 
tohydrodynamical simulations' made under the European 
Heads of Research Councils and European Science Founda- 
tion EURYI (European Young Investigator) Awards scheme, 
was supported by funds from the participating organizations 
of EURYI and the EC Sixth Framework Programme. 



6 CONCLUSIONS 

We have carried out radiation hydrodynamical simulations 
to investigate the evolution of massive self-gravitating discs 
(Afdisc/Af* = 0.1). We consider discs with opacities rang- 
ing from 0.01 X to lOx the interstellar Rosseland mean val- 
ues. We also consider the effects of changing the initial and 
boundary temperatures of the discs as well as simulating 
different disc sizes (with outer disc radii, -Rout = 25 and 300 
AU). 

We find that the disc opacity is very important in de- 
termining whether a disc is likely to fragment. In particular, 
we find that fragmentation is promoted in discs with opacity 
values lower than interstellar Rosseland mean values since 
this allows radiation to leave the disc quickly. Low opaci- 
ties may exist in low metallicity discs or discs with larger 
grain sizes. This is a particularly important and timely re- 



sult given the recent discoveries of wide orbit planets ( Kalas 
Marois et al. 2008) and the future emphasis 



et al. 2008 



for surveys of planets on such wide orbits. We show that 
it is possible for fragmentation to occur in gravitationally 
unstable discs even at radii where the innermost planet of 
the HR 8799 system is located (R > 24 AU). Furthermore, 
HR 8799 is known to be a metal-poor, A Bootis star with 
metallicity [M/H] = -0.47 ( |Gray fc Kaye|1999| ) so it is rea- 
sonable to assume that its disc was similarly metal-poor. We 
have shown that such a scenario favours fragmentation and 
therefore, our results indicate that all three planets of the 
HR 8799 system may well have formed via gravitational in- 
stability. Though a hybrid core accretion and gravitational 
instability scenario for planet formation may also be a pos- 
sibility for this system, our calculations show that such a 
hybrid scenario may not be necessary. 

We find that the presence of a thermal blanket as a re- 
sult of the stellar irradiation inhibits fragmentation since 
the discs are only able to cool to the boundary temperature. 
However, we also show that under certain circumstances, 
fragmentation may occur. Our results demonstrate that for 
fragmentation, weak irradiation is required such that the 
boundary temperature and hence Toomre stability parame- 
ter is low, in addition to low enough opacities (even in large, 
cool discs) since this allows more efficient cooling so that the 
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